Separatrices and basins of stability from time series data 
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An approach is presented for identifying separatrices in phase space generated from noisy time 
series data sets representative of measured experimental data. These separatrices are identified as 
ridges in the phase space distribution of finite-time Lyapunov exponents, i.e., Lagrangian coherent 
structures (LCS). As opposed to previous approaches, the LCS is identified using only trajectories 
since no analytical or data-defined vector field is available. The method is applied to a biological 
simulation in which the separatrix reveals a basin of stability. These results suggest that the method 
will be a fruitful approach to time series analysis, particularly in cases where a limited number of 
trajectories are available as might be encountered in experiments. 
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Increasingly, dynamical systems of interest are defined 
not by analytical models, but by data from experiments 
or large-scale simulations, such as musculoskeletal biome- 
chanics [U [2j [3] or geophysical fluid dynamics [U [5j EJ [7] . 

In many cases, researchers want to ascertain if deter- 
ministic chaos is present [2 , 8 . This can be acheived by 
determining characteristic exponents that describe the 
sensitivity of the solution to initially close starting con- 
ditions. One popular technique is to estimate the (max- 
imum) Lyapunov exponent averaged over the sampled 
portion of phase space [TJ3 EH H3 E2 H3 . This method 
is well suited for analysis of time series data from exper- 
iments. 

However, there is more information contained within 
the maximum Lyapunov exponents. Instead of averag- 
ing the exponents over phase space to obtain a single 
scalar value, one can generate a maximum Lyapunov ex- 
ponent field by considering how phase space expansion 
is distributed. From this field, one can determine sep- 
aratrices, co-dimension one boundaries in phase space 
separating qualitatively different kinds of motion. 

An example from biomechanics will demonstrate our 
point. In biomechanics, a separatrix exists between 
standing and falling. Standing with postural sway is a 
distinctly different type of motion than falling. During 
standing, the body remains in the vicinity of an equilib- 
rium position and may be characterized as dynamically 
stable. Compare this motion to falling where the body 
rapidly diverges from the equilibrium position at an in- 
creasing velocity. In falling, the body behaves unstably 
with respect to the upright vertical position. 

If we allow an experimental subject to take a step dur- 
ing fall recovery, another boundary will develop. Now 
three states exist: standing, recovering from a fall with 
one step, and falling. Each type of motion is divided from 
the other by a separatrix. Extending this theory, a phase 
space diagram with multiple fronts may be generated. 

Similar to standing postural sway is the problem of 
maintaining torso stability. Torso stability is necessary 
to avoid large deformations in the lumbar spine that may 
result in low back injury and pain. In this Rapid Commu- 



nication, a method will be developed to locate the basin 
of stability that separates stable, injury free torso sway 
from unstable, potentially injurious motion. 

The computation of finite-time Lyapunov exponents 
from experimental data has been used before in muscu- 
loskeletal biomechanics, particularly to quantify local dy- 
namic stability during locomotion [1 . We demonstrate 
that the phase space distribution of finite-time Lyapunov 
exponents can also provide information. This informa- 
tion includes the boundaries, or separatrices, between 
qualitatively different states, as described above for the 
biomechanics example. 

The separatrices are found as Lagrangian coherent 
structures (LCS), borrowing a term from the fluid me- 
chanics [15] , and are defined as the ridges of the finite- 
time Lyapunov exponent field [6 j. These structures indi- 
cate the location of the separatrix demarking the bound- 
ary between qualitatively different kinds of motion. LCS 
have previously been used to analyze dynamical systems 
defined by fluid flow fields from data [16] [17] and analyt- 
ical biochemical models [18] but have not been applied 
to biomechanical data or any time series data obtained 
from experiments. 

The goal of this Rapid Communication is to demon- 
strate that LCS can be used to analyze experimental 
time series data where only individual trajectories are 
available. LCS are robust with respect to noise [16] [T9] 
making them even more attractive for use in experimen- 
tal data analysis where noise sensitivity is an important 
issue [201 ED 122] • 

This Rapid Communication is organized as follows: 
first, we discuss the method for estimating the maximum 
finite-time Lyapunov exponent (FTLE) associated with 
a point in phase space. Second, we demonstrate how 
the phase space distribution of maximum FTLE can be 
obtained using only trajectories. Third, using a model 
of human postural control, we demonstrate how ridges 
can be found in the phase space distribution of the max- 
imum FTLE. These ridges are separatrices partitioning 
the phase space into regions of distinct behavior. 

In the following discussion, we borrow ideas developed 
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in fluid mechanics [15] US H3 US [24] and with some 
modifications apply them to a new context, a situation 
where only trajectories are available and not the vector 
field itself. 

Sensitivity analysis and finite-time Lyapunov 
exponents. — Suppose we are given a reference tra- 
jectory x(t) going from xq at time to to x\ at time t\. 
We assume the trajectory evolves under the dynamical 
equations of a time-independent (autonomous) system 



x = f(x), xe 



(1) 



This equation describes a flow field or vector field. The 
sensitivity of the reference trajectory is discussed below. 

Let trajectories of the system ([I) with x(to) = xo be 
denoted by </>(t,to). In other words, 0(t, to) : x(to) i— > 
x(t) denotes the flow map of the dynamical system ([!]), 
mapping particles from their initial location at time to 
to their location at time t. For our purposes, we will 
denote the flow map as (j)(t,to]Xo) or simply <j)(t]Xo) so 
the dependence on the initial condition x(to) — xo is 
made clear. 

Consider a second trajectory that starts slightly away 
from the reference trajectory x(t), i.e., starts from the 
perturbed initial vector xo + Sxo at time to. As the tra- 
jectories evolve, the vector displacement (or perturbation 
vector) 



Sx(t) = </>(t; xo + Sxq) — </>(t; xo) 



(2) 



will also evolve. For our purposes, the "second trajec- 
tory" might be the result of another experimental trial 
or another portion of the same trajectory separated in 
time by a minimum threshold value. We discuss this fur- 
ther below. 

The linear relationship between small initial perturba- 
tions and perturbations at some time t is 



Sx(t) = $(t, t )feo- 



(3) 



d<t>(t;xo) 
dx 



is the state transition matrix 



where $(Mo) 

(also known as the fundamental matrix). The state tran- 
sition matrix can be viewed as a deformation gradient. 
If an (infinitesimal) n-dimensional spherical blob of par- 
ticles is placed about the reference trajectory, then af- 
ter a duration T = t — to, the blob will have expanded 
in some directions and compressed in others to form an 
n-dimensional ellipsoid. The matrix 3>(t, to) contains in- 
formation about this expansion and contraction as well 
as the rotation of the initial blob of particles, due to the 
locally deforming nature of the flow. 

Suppose there exists a state transition matrix over 
some interval, 3>(t, to)- The size of the final perturba- 
tion at time t is given by 



\\5x(t)\\ 2 = 5x* Mt,to)*$(t,t )}6x , 



(4) 



where || • || is the vector norm on R n , A* denotes the 
transpose of the matrix A, and the perturbations are 
considered as column vectors. The symmetric matrix 



is the finite-time right Cauchy-Green deformation tensor 
[6]. The matrix C is a rotation-independent measure 
of deformation; it gives the square of the local change 
in distances due to deformation [25j [26]. Since C is a 
symmetric, positive definite matrix, it has n real, positive 
eigenvalues [271128] . 

One can associate with point xo a maximum finite-time 
Lyapunov exponent, given by 



dxixo) 



c(C), 



(6) 



where T = t — to is the finite duration over which ex- 
pansion is measured and A max (C) is the maximum eigen- 
value of C with the corresponding (normalized) eigenvec- 
tor ei(to). In other words, if Sxo is along ei(to) at time 
to, then maximum stretching occurs over the time T and 
the length of the perturbation vector becomes 



\\5x{ 



= e <xiOo)T 



(7) 



C = $(Mo)*$(Mo), 



(5) 



where t = t + T| 

Computing the maximum FTLE field when only indi- 
vidual trajectories are available. — If we assume that the 
direction of maximum expansion dominates the dynam- 
ics of perturbations in arbitrary directions [13], we can 
approximate the maximum FTLE field. Under this as- 
sumption, we assume that (7|) holds for all perturbation 
vectors, regardless of initial phase space direction. This 
will result in a conservative approximation to the actual 
value when the arbitrary vector direction is not aligned 
with the direction of maximum expansion. As our con- 
cern is now focused on the maximum FTLE, we will refer 
to the maximum FTLE as simply the FTLE hereafter. 

The FTLE is estimated as the rate of separation of 
neighboring trajectories. In order to understand how this 
is determined from experimental data, consider the ref- 
erence trajectory shown in Figure [I] With a reference 
point established, a target location (pi) is identified that 
is a perturbation distance Sq from the reference point in 
phase space. The data point closest to pi on another 
trajectory is then found, n\. The other trajectory can 
be from either a different experimental trial or another 
portion of the same trial separated by a sufficient amount 
of time to avoid a correlation with the reference point. 
This process is repeated for other directions. Using this 
method, 2n neighbors are found for each reference point 
corresponding to positive and negative directions of each 
dimension of state space. If multiple nearest neighbors 
are considered, perturbations are sampled in multiple 
phase space directions which increases the likelihood of a 
separation lying in the direction of maximum expansion. 

Each point in the data set is sequentially evaluated 
by considering it to be a reference point. The FTLE is 
calculated for each pair, but unlike previous analyses that 
averaged the FTLE over time and phase space [2j [3j [TTJ 
[T3] , the exponent value over a finite time T is associated 
with a phase space location midway between the pair 
of points. In this way, the phase space distribution of 
the FTLE is developed. The FTLE field is generated 
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FIG. 1: Estimating the maximum FTLE by averaging over the 
growth of perturbation vectors in multiple phase space directions. 
We make the assumption that the maximum FTLE dominates the 
evolution of the perturbation vectors. 



by placing an n-dimensional grid over the phase space 
and using the distribution of FTLE to generate an n- 
dimensional surface. At each phase space location, x, 
the height of the surface is the scalar value cri(x) from 
equation^. 

The choice of Sq influences the FTLE field. The pertur- 
bation distance 5q is a coarse-graining parameter selected 
to be large enough to overcome system noise and small 
enough to reveal local features of the FTLE field. 

Figure [2] is used to explain why ridges develop in the 
FTLE field at the separatrices. In this schematic rep- 
resentation, divergence of two points (a and b) within 
the stable region and two points (c and d) in the unsta- 
ble region are compared. As the points a and b evolve 
in time, their trajectories only slightly diverge resulting 
in a small FTLE. This value forms a data point on the 
FTLE field midway between a and b. Similarly, points 
c and d also diverge slightly resulting in a small FTLE 
and another data point for generating the FTLE field. 
However, b and c on opposites sides of the separatrix 
diverge greatly even over short times. This results in a 
large FTLE that forms a data point on the FTLE field 
at the approximate location of the separatrix. As more 
data points are entered a "volcano shaped" ridge forms 
revealing the separatrix. 

Numerically simulated experiments with noise. — We 
have implemented the above approach for a model of hu- 
man postural control: the inverted pendulum with lim- 
ited gain proportional control. A simulation with 20 in- 
dependent trials was generated from the model (Figure 
[3|. In the simulation, the system begins at equilibrium 
and is perturbed with Gaussian random force perturba- 
tions. These forces generate movement that is attenuated 
by a controller. As a result, the system is able to main- 
tain stability for a period of time before an unrecoverable 
sequence of perturbations causes the system to become 
unstable. Therefore, an important aspect of this simu- 
lation is that it spans the space of possible states. This 
is important because it is not possible to experimentally 
determine the edge of the basin of stability unless data 
is available on both sides of the separatrix. 




FIG. 2: The divergence of two points of opposite sides of a separa- 
trix is larger than the divergence of points on the same side. This 
generates a ridge in the FTLE field at the separatrix. 



Phase space reconstruction is considered a preprocess- 
ing step in our method. For the biomechanical system 
used to demonstrate the method, we do not use time- 
delay reconstruction. Since the system under study is a 
mechanical system, we take the measured coordinate (q) 
and numerically construct the time derivative (g), taking 
the 2D space of x = (g, q) as the reconstructed phase 
space. However, we note that the method of using an 
FTLE field to find separatrices is not tied to any partic- 
ular means of phase space reconstruction. 

The resulting FTLE field and separatrix for an evolu- 
tion time of T = 1.2 seconds are shown (Fig. [4|. This 
value of T approximately corresponds to the characteris- 
tic timescale of evolve in the unstable region. 

Two ridges in the distribution of maximum FTLE are 
discernible. As shown elsewhere [6], ridges in the FTLE 
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FIG. 3: The time-series data analyzed came from several simulated 
experimental trials. 
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FIG. 4: A clearly observable separatrix forms the boundary of the 
basin of stability around an equilibrium region. 

field reveal the location of partial transport barriers or 
separatrices in the phase space, providing boundaries 
which partition the space into regions of different be- 
haviors. In this case, the separatrix forms a boundary 
between the region of stable postural sway (around the 
origin) and falling motion (beyond the boundary). The 
characteristics of the stable region and its boundary loca- 
tion in this biological example depend on the accuracy of 
a number of neurological sensory systems, the feedback 
gain associated with core muscle strength, and the time 
delay of the postural control system. 

Although 20 trials were evaluated for this simulation, 
a separate simulation indicates that LCS structures may 
be identified with as little as two time series using this 



method. However, smaller trial numbers will tend to be 
more affected by random noise than larger sets because 
the mean value of random noise approaches zero as the 
number of trials increases. 

Conclusions and future directions. — We have shown in 
this Rapid Communication that even without an analyt- 
ically defined or data defined vector field, we are able 
to identify separatrices using the LCS method. A dy- 
namic model of human postural control driven by noise 
was used to generate a time series similar to experimen- 
tal data. When we applied the LCS method to the sim- 
ulated data, a separatrix was revealed which formed the 
boundary of a basin of stability around an equilibrium 
location. As a result, we believe this method provides a 
fruitful approach for extracting information from noisy 
experimental data on boundaries between qualitatively 
different kinds of motion. 

We note that the phase space averaged FTLE for any 
time T can be obtained by computing the average of the 
FTLE field over the sampled region of phase space. This 
provides the link between the current method and previ- 
ous methods for finding an averaged Lyapunov exponent. 

In forthcoming work, we will demonstrate the method 
on higher dimensional data and actual experimental data. 
We will also consider another method for constructing 
an FTLE field, one which estimates the state transition 
matrix $(t,to). We believe this may be a more accurate 
method for obtaining the FTLE field over short evolution 
times, T. 
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